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ABSTRACT 

A two-dimensional (2D) shock-rest-frame model for particle simulations is developed. Then full kinetic 
dynamics of a perpendicular collisionless shock is examined by means of a 2D full particle simulation. We 
found that in the 2D simulation there are fewer nonthermal electrons due to surfing acceleration which was 
seen in the previous ID simulations of a high Mach number perpendicular shock in a low -beta and weakly 
magnetized plasma. This is because the particle motion along the ambient magnetic field disturbs the formation 
of coherent electrostatic solitary structures which is necessary for electron surfing acceleration. 
Subject headings: acceleration of particles — methods: numerical — shock waves — supernova remnants 



1. INTRODUCTION 

Observations of shell-type supernova remnants (SNRs) pro- 
vide us evidence for electron acceleration at SNR shocks, and 
in some c ases the maximum energy of the electrons ex ceeds 
~10 TeV (iKoyama et alJl99llBamba et al.l2003Ll2005l) . The 
most plausible mechanism of the acceleration is the diffu- 
sive shock acceleration (DSA), which explains t he broadband 
power-law spectrum with an index around 2 (iDruryl IT983b 
iBlandford & EichlerlfT987h . Before the DSA phase in which 
electrons cross the shock front many times, they have to be 
pre-accelerated by an unknown injection mechanism. This 
"injection problem" is still unresolved. 

One of the possible inje ction mechanisms i s the 
shock surfing acceleratio n dShimada & Hoshinol 120001 : 
iHoshino & Shimadal 120021) . IShimada & Hoshinol j2000l) 
performed a one-dimensional (ID) electromagnetic full par- 
ticle simulation and found the formation of large-amplitude 
electrostatic solitary structures during the cyclic reformation 
of a high Mach number perpendicular shock in a low-beta and 
weakly magnetized plasma. The coherent solitary structures 
trap electrons in their electrostatic potential well, resulting 
in significant surfing acceleration of electrons. The electron 
surfing has been regarded as an efficient direct acceleration 
mecha nism at SNR shocks and studied by many authors 
(e.g.. iMcClements et all 120011: IHoshino & Shimada 120021: 
ISchmitz et al]|2002allbt lAmano & Hoshinoll2007l) . However, 
these previous works are based on ID simulations, in which 
configuration of shock magnetic fields cannot be modified. 
Since charged particles can easily drift along magnetic fields, 
it is important to take into account the modification of shock 
magnetic fields and the free motion of particles along shock 
magnetic fields. Hence, multi-dimensional simulations are 
necessary. 

In the present study we perform a 2D full particle simula- 
tion of a high Mach number and low -beta perpendicular colli- 
sionless shock. The preliminary result shows that the electron 
surfing acceleration is suppressed in a 2D system. 



2. FULL PARTICLE SIMULATION 



1 Solar-Terrestrial Environment Laboratory, Nagoya University, Nagoya, 
464-8601, Japan; umeda@stelab.nagoya-u.ac.jp 

2 Department of Physical Science, Hiroshima University, Higashi- 
Hiroshima 739-8526, Japan; ryo@theo.phys.sci.hiroshima-u.ac.jp 



Generally speaking, it is difficult to perform multi- 
dimensional simulations of collisionless shocks even with 
present-day supercomputers. In the previous simulations of 
high Mach number s ho cks (e.g., Shimada & Hoshinol 120001: 
SchmitzetalJ l2002allbl : iLeeetalJ I2004F lAmano & Hoshinol 
2007), a collisionless shock is excited with the "injection 
method" (also called the reflection, wall, or piston method), 
in which a plasma is injected from one side of the simula- 
tion domain and is reflected back when it reaches the other 
side. A collisionless shock is excited by interaction between 
injected and reflected plasmas. Thus, the simulation domain 
is taken in the downstream rest frame, and the excited shock 
wave propagates upstream. It is necessary to take a very long 
simulation domain in the propagation direction of the shock 
wave in order to study long-time evolution of moving shock 
waves. 

Recently, ID full particle si mulations have been performed 
in an almost shock rest fr ame (Muschi eitti & Lembegdl2006T : 
lUmeda & YamazakH 120061) . In these simulations, a col- 
lisionless shock is excited with the "relaxatio n method," 
which was fi rst u sed in hybrid code simu la tions dLeroy et all 
I1981L 119821) . iMuschietti & Lembegd d2006l) used the 
Darwin approximation in which the displacement current 
due to the transve r se ele ctric field is neglected, while 
lUmeda &~Y amazaki (2006) solved the full set of Maxwell's 
equations to take into account all the plasma waves. Although 
the shock-rest-frame model is useful to perform simulations 
of shocks with less computer resources, there have been few 
multi-dimensional shock-rest-frame models. Here, we suc- 
cessfully extended the pre vious ID shock-rest-frame model 
(lUmeda & Yamazakill2006l) to two dimensions; 

We use a 2D electromagnetic particle code ( Umeda 2004), 
in which the full set of Maxwell's equations and the relativis- 
tic equation of motion for individual electrons and ions are 
solved in a self-consistent manner. The continuity equation 
for charge is also solved to compute the exact current densit y 
given by the motion of charged particles (lUmeda et alj|2003l) . 

The initial state consists of two uniform regions separated 
by a discontinuity. In the upstream region that is taken in 
the left-hand side of the simulation domain, electrons and 
ions are distributed uniformly in space and are given random 
velocities (v x ,v y ,v z ) to approximate shifted Maxwellian mo- 
mentum distributions with the drift velocity u x \ , number den- 
sity n\ = eom e oj^ el /e 2 , and isotropic temperatures T e \ = m e vf el 
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and T,i = fn(vf a , where to, e, w p , and v t are the mass, charge, 
plasma frequency, and thermal velocity, respectively. Sub- 
scripts "1" and "2" denote "upstream" and "downstream," re- 
spectively. The upstream magnetic field Byoi = —m e uj ce i/e is 
also assumed to be uniform, where ui c is the cyclotron fre- 
quency (with sign included). The downstream region taken 
in the right-hand side of the simulation domain is prepared 
similarly with the drift velocity u X 2, density m, isotropic tem- 
peratures T e 2 and 7^, and magnetic field B v o2- 

We take the simulation domain in the x-y plane and as- 
sume a perpendicular shock (i.e., B x o = 0). Since the ambi- 
ent magnetic field is taken in the y-direction, free motion of 
particles along the ambient magnetic field is taken into ac- 
count. As a motional electric field, a uniform external electric 
field E z o = —UxiByQi = —u&Byoz is applied in both upstream and 
downstream regions, so that both electrons and ions drift in 
the x-direction. At the left boundary of the simulation domain 
in the x-direction, we inject plasmas with the same quantities 
as those in the upstream region, while plasmas with the same 
quantities as those in the downstream region are also injected 
from the right boundary in the x-direction. We adopted ab- 
sorbing boundaries to suppress nonphysical reflection of elec- 
tromagnetic w aves at both ends o f the simulation domain in 
the x-direcfion (lUmeda et al.l200 lb , while the periodic bound- 
aries are imposed in the y-direction. 

The initial downstream quantities are given by solving the 
shock jump conditions for a magnetized tw o-fluid isotropi c 
plasma consisting of electrons and ions (e.g.. lHudsonF l970). 
In order to determine a unique initial downstream state, we 
need given upstream quantities u x i, uj pe \, Wcei, VteU an d 
Tn/T e \ and another parameter. We assume a low-beta and 
weakly magnetized plasma such that (3 e \ = (3n = 0.125 and 
uj ce i/w pe \ = -0.05 in the upstream region. We also use a 
reduced ion-to-electron mass ratio m;/m e = 25 for compu- 
tational efficiency. The light speed c/v, e \ = 80.0 and the 
bulk flow velocity of the upstream plasma u x \/v te \ = 8.0 are 
also assumed. Then, the Alfve n Mach number is calcu- 
lated as Ma = (u x \ / c)\lo pe \ / w ce \\^J mi / m e = 10.0. The ion-to- 
electron temperature ratio in the upstream region is given as 
Tn/Tei = 1.0. Note that these parameters are almost the same 
as the previous full particle simulations of a perpendicular 
shock ( Shimada & Hoshino 2000; Hosh ino & Shimadall2002t 
ISchmitz et all l2002allbt iLee et all 120041) . In this study, the 
downstream ion-to-electron temperature ratio TnjT e 2 = 8.0 
is also assumed as an initial parameter to obtain the unique 
downstream quantities ui pe 2/ui pe \ = 1.953, u> C e2/wpei = 0.191, 
ua/vtei =2.097, and v te 2/v, e \ =5.715. 

We used N x x N y = 4096 x 64 cells for the upstream re- 
gion and N x x N y = 2048 x 64 cells for the downstream re- 
gion. The grid spacing and time step of the present simula- 
tion are Ax/\ el = 1.0 and ui pel At = 0.005. We used 16 pairs 
of electrons and ions per cell in the upstream region and 64 
pairs of electrons and ions per cell in the downstream region 
at the initial state. One might think that the downstream sys- 
tem size is too short because the downstream does have an 
affect on the transition region. However, we have confirmed 
by changing the downstream system size that there is not any 
critical change in the simulation result when the downstream 
system size is longer than t he typical ion gyroradius u x \ / uj c a 
(Umeda & Yamazaki 2006). This implies that the results are 
most appropriate to the transition region, and that the transi- 
tion region, is approximately modeled by the present shock- 
rest-frame model despite the lack of an extensive downstream 
region. Thus, the processes seen in the present simulation 
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FIG . 1 . — Transverse magnetic field By as a function of position x and time 
t on the .v-axis (y = 0). The position and time are normalized by p, = u x i/w c n 
and 1 /w c n , respectively. The magnitude is normalized by the initial upstream 
magnetic field B,,oi- 



should be robust. 



3. RESULT 



Figure 1 shows the transverse magnetic field B v as a func- 
tion of position x and time t on the x-axis (y = 0). In the 
present shock-rest-frame model, a shock wave is excited by 
the relaxation of the two plasmas with different quantities. 
Figure 1 shows that the shock front appears and disappears 
at a timescale of the downstream ion gyro-period, which cor- 
responds to the cyclic reformation of high Mach number per- 
pendicular shocks. Since the initial state is given by the shock 
jump conditions for a two-fluid plasma consisting of electrons 
and ions, the excited shock is "almost" at rest in the simula- 
tion domain. The shock front shifts upstream at a very slow 
velocity of ~ 0.06m v i, because dynamics of reflected ions is 
not taken into account in the two-fluid approximation. We 
have also performed a ID simulation with the same parameter 
as the 2D simulation and confirmed that the period of cyclic 
reformation in both simulations is the same. 

Figure 2 shows the shock-normal electric field E x as a func- 
tion of position x and time t on the x-axis (y = 0). There is a 
good correlation between Figures 2 and 1, implying that the 
behavior of a shock-normal electric field is controlled by the 
shock reformation process. Figure 2 shows a bipolar signature 
of shock-normal electric field turning from negative to posi- 
tive in the spatial scale of ion gyroradius from x/p, = -0.2 
to 0.2. The ion-scale electric field corresponds to the shock 
potential which reflects ions upstream. In the shock foot re- 
gion there exist sinusoidal wave structures, which are excited 
by the current-driven instability due to reflected ions. These 
features are also seen in the ID simulation (not shown here). 



2D PIC Simulation of a Perpendicular Shock 
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FIG. 2. — Shock-normal electric field Ex as a function of position x and time 
t on the .Y-axis (v = 0). The position and time are normalized by p, = u x \ /LO cji 
and \/u c n, respectively. The magnitude is normalized by the motional elec- 
tric field E z q. 

However, there are two differences between the ID and 2D 
results. The first is that the amplitude of electric field in 
the 2D simulation is much smaller than that in the ID sim- 
ulation. The second is that the ID simulation shows more 
complex wave structures with several large-amplitude pulses. 
The typical amplitude of coherent solitary structures in the 
previous ID simulations was E x / E 7 q ~ 20 and the maximum 
amplitude exceeded E x /E 7 q ~ 40 (IShimada & Hoshinoll2000t 
iHoshino & Shim ada 2002). In the present 2D simulation, 
on the other hand, the typical amplitude of a shock-normal 
electric field is E x /E z q ~ 5 and the maximum amplitude is 
E x /E z o ~ 10. This result suggests that the saturation level of 
the current-driven instability in the present 2D simulation is 
much lower than that in the previous ID simulation. 

Figure 3 shows snapshots of position-velocity phase-space 
distribution functions of electrons and ions at u) c nt = 4.91 and 
the corresponding spatial profiles of E x and B x components. 
At this time, there is a strong ion beam reflected upstream. 
Note that the x—v x phase-space distribution functions are re- 
duced by integrating over y, v y , and v z , and the x - v y phase- 
space distribution function is reduced by integrating over y, 
v x , and v.. 

Strong heating of electrons is found from the phase-space 
plots in Fig. 3 in the shock-normal direction in the shock foot 
region (xj p-, = -0.2-0) and at the shock overshoot (x/p,- = 0. 1- 
0.3). The shock-normal electric field E x is strongly enhanced 
by the current-driven instability due to reflected ions in the 
shock foot region. An electron phase-space vortex is also 
identified at the leading edge of reflected ions (x/pi ~ —0.2). 
These results are consistent with the previous ID simulations. 

The beam-plasma interaction in the shock-normal direction 
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FIG. 3. — Position- velocity phase-space distribution functions of electrons 
and ions in x — v xe , x — v ye , and x— v x i phase spaces, at u) c nt = 4.91, and the 
corresponding spatial profiles of E x and B x components. Note that the x—v x 
phase-space distribution functions are reduced by integrating over y, v v , and 
v z , and the x— Vy phase-space distribution function is reduced by integrating 
over y, v x , and v z . 

results in strong temperature anisotropy T e ±/T e u > 10. We 
found enhancement of the shock-normal magnetic compo- 
nent B x especially in the electron phase-space vortex and in 
the shock foot region, which corresponds to electromagnetic 
whistler mode waves due to electron temperature anisotropy. 
However, the amplitude of the whistler mode is not so strong 
that the structure of the transverse magnetic field B v is not 
modified, keeping an almost ID structure. 

In the electron x— v v phase-space plot, the strong heating of 
electrons is found in the electron phase-space vortex, in the 
shock foot region, and at the shock overshoot. These elec- 
trons are thermalized along the ambient magnetic field by the 
whistler mode wave. The electron temperature anisotropy, en- 
suing excitation of electromagnetic whistler mode waves, and 
electron heating along the ambient magnetic field would be 
common features in perpendicular shocks. However, these 
were not seen in the previous ID simulation. 

Figure 4 shows energy distribution functions of electrons 
in the downstream region for ID and 2D simulations. The 
solid line shows the 2D simulation result, and the dashed 
line shows the ID simulation result. It is known that there 
appear nonthermal electrons abo ve the upstream bulk en- 
ergy of ions in ID simul ations ([Shimada & Hoshino 120001 : 
IHoshino & Shim ada 2002). On the other hand, we found less 
nonthermal electrons in the present 2D simulation. 

The previous ID studies reported electron surfing accelera- 
tion with coh erent electrostatic solitary structures in the tran- 
sition regio n ( Shimada & Hoshino 2000; Hoshino & Shimada 
I2002t ISchmitz et alj|2002allbl : lUmeda & Yamazakill2006l) . In 
the present 2D simulation, however, we could not confirm the 
existence of solitary structures in the shock foot region. It is 
well known that coherent electrostatic solitary structures gen- 



4 



Umeda, Yamao & Yamazaki 



Electron energy spectin 




FIG. 4. — Energy distribution functions of electrons in the downstream 
region for ID and 2D simulations. The solid line shows the 2D simulation 
result and the dashed line shows the ID simulation result. The dotted line 
indicates Maxwellian distribution. The energy is normalized by the upstream 
bulk energy of ions. 

erated by nonlinear beam-plasma interactions need a strong 
guiding magnetic field to maintain a n electrostatic equilib - 
rium state, the so-called BGK mode dBernstein et al.|[l957l) . 
In the present simulation model, electrostatic fields appear in 
the x-direction, while a weak ambient magnetic field is taken 
in the y-direction. Since there is no guiding magnetic field in 
the x-direction, electrons can freely move along the ambient 
magnetic field, by which electrostatic waves diffuse. As a re- 
sult, the saturation level of electrostatic waves becomes much 
lower than that with a strong guiding magnetic field (e.g., 
Morse & Nielsonlll969t iMuschietti et al..ll2000t lUmeda et alj 
2006), and the excited electrostatic waves do not develop into 
coherent solitary structures. 

In addition, waves can also propagate along the ambient 
magnetic field in the present 2D model. Electromagnetic 
whistler mode waves are excited by temperature anisotropy, 
which is caused by electron heating in the shock-normal di- 
rection due to current-driven instability. The excited whistler 
mode waves propagate along the magnetic field and thermal- 
ize electrons in the direction parallel to the ambient mag- 
netic field. As a result, we found strong heating of electrons 
along the magnetic field instead of electron surfing acceler- 
ation across the magnetic field. We did not find any further 
acceleration by the whistler mode waves. 



4. SUMMARY AND DISCUSSION 

We developed a 2D shock-rest-frame model for full parti- 
cle simulations. Then a 2D simulation of a low-beta and high 
Mach number perpendicular shock was performed, which 
shows that electrons are thermalized in the direction paral- 
lel to the ambient magnetic field in contrast to the previous 
ID simulations. The heating of electrons along the ambient 
magnetic field is due to electromagnetic whistler mode waves 
excited by temperature anisotropy. In addition, the current- 
driven instability due to reflected ions does not lead to for- 
mation of coherent electrostatic solitary structures because of 
the free motion of electrons along the ambient magnetic field. 
As a result of these consequences, the generation of nonther- 
mal electrons via electron surfing acceleration due to electro- 
static solitary structures is suppressed in the 2D system. This 
process is independent of shock parameters such as uj pe /u) ce , 
(3, and Ma- Even with the realistic parameters, there would 
be less surfing acceleration in purely perpendicular shocks, 
because there is no guiding magnetic field to stabilize coher- 
ent electrostatic structures. However, the surfing acceleration 
might be effective in an oblique shock in which a guiding 
magnetic field could allow the formation of large-amplitude 
electrostatic struc tures. 

It is noted that lOhira & Takaharal d2007l) have also shown 
that electron surfing acceleration becomes absent by a sim- 
plified 2D simulation with different shock geometry. That is, 
their simulation domain has been taken in a plane perpendic- 
ular to the ambient magnetic field (which corresponds to the 
x-z plane in the present model). Thus, electrostatic solitary 
structures generated by the current-driven instability lose co- 
herency by a quite different mechanism. However, both 2D 
simulations with different axes of coordinate have indicated 
that the electron surfing acceleration is suppressed more than 
in ID simulations. Hence, another unknown injection mecha- 
nism for electron acceleration might be necessary. 
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